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Abstract. We provide a real algebraic symbolic-numeric algorithm for computing the real va- 
riety Vr(7) of an ideal / C K[x], assuming Va(T') is finite (while Vc(I) could be infinite). Our 
approach uses sets of linear functionals on M[x], vanishing on a given set of polynomials generat- 
ing / and their prolongations up to a given degree, as well as on polynomials of the real radical 
ideal \/7 obtained from the kernel of a suitably defined moment matrix assumed to be positive 
semidefinite and of maximum rank. We formulate a condition on the dimensions of projections 
of these sets of linear functionals, which serves as stopping criterion for our algorithm; this new 
criterion is satisfied earlier than the previously used stopping criterion based on a rank condition 
for moment matrices. This algorithm is based on standard numerical linear algebra routines and 
semidefinite optimization and combines techniques from previous work of the authors together 
with an existing algorithm for the complex variety. 



1. Introduction 

Polynomial equations play a crucial role in mathematics and are widely used in an emerging 
number of modern applications. Recent years have witnessed a new trend in algebraic geometry and 
polynomial system, namely numerical polynomial algebra [26] or numerical algebraic geometry [25]. 
Algorithms in this field deal with the problem of (approximately) computing objects of interest in 
the classical area of algebraic geometry with a focus on polynomial root finding. 

There is a broad literature for the problem of computing complex roots, that deals with numerical 
and symbolic algorithms, ranging from numerical continuation methods as in e.g. Verschelde [28] 
to exact methods as in e.g. Rouillier [23], or more general Grobner or border bases methods; see 
e.g. the monograph [9] and the references therein. 

In many practical applications, one is only interested in the real solutions of a system of poly- 
nomial equations, possibly satisfying additional polynomial inequality constraints. An obvious 
approach for finding all real roots of a system of polynomial equations is to first compute all com- 
plex solutions, i.e., the algebraic variety Vc(I) of the associated ideal / C M[x], and then to sort 
the real variety Vr(I) = R n D Vc(-0 from Vc(I) afterwards. However, in many practical instances, 
the number of real roots is considerably smaller than the total number of roots and, in some cases, 
it is finite while |Vc(-0| = oo. 

The literature about algorithms tailored to the problem of real solving systems of polynomial 
equations is by far not as broad as for the problem of computing complex roots. Often local Newton 
type methods or subdivision methods based on Descartes rule of sign, on Sturm-Habicht sequences 
or on Hcrmite quadratic forms are used; see e.g. [1, 20, 22] for a discussion. In [12] we gave an 
algorithm for finding Vr(I) (assumed to be finite), and a semidefinite characterization as well as 
a border (or Grobner) basis of the real radical ideal \fl : by using linear algebra combined with 
semidefinite programming (SDP) techniques. We exploited the fact that all information needed to 
compute the above objects is contained in the so-called moment matrix (whose entries depend on 
the polynomials generating the ideal /) and the geometry behind when this matrix is required to be 
positive semidefinite with maximum rank. We use the name (real-root) moment-matrix algorithm 
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for the algorithm proposed in [12]. This algorithm was later extended to the computation of all 
complex roots in [13]. A feature of the real-root moment-matrix algorithm is that it requires 
to solve a sequence of SDP problems involving matrices of increasing size until a certain rank 
condition is satisfied. Solving the SDP problem is the computationally most demanding task in 
the algorithm. It is thus important to be able to terminate the algorithm at an as early as possible 
stage so that the size of the matrices does not grow too much. This is the motivation for the 
present paper where we present a new stopping condition, which is satisfied at least as early as the 
rank condition of [12] (and often earlier on examples). This leads to a new algorithm which we 
name (real-root) prolongation-projection algorithm since its stopping condition involves computing 
the dimensions of projections of certain sets of linear functionals on spaces of polynomials. This 
new algorithm arises by incorporating several ideas of [12, 13] into an existing symbolic-numeric 
solver dedicated to compute Vc{I) (as described e.g. in [32]). A detailed description will be given 
in Section 5 but, in order to ease comparison with the moment-matrix method of [12], we now give 
a brief sketch of both methods. 

Sketch of the real-root moment-matrix and prolongation-projection algorithms. While 
methods based on Grobner bases work with the (primal) ring of polynomials R[x], its ideals and 
their associated quotient spaces, we follow a dual approach here. The algorithms proposed in [12] 
and in this work manipulate specific subspaces of (R[x])*, the space of linear forms dual to the 
ring of multivariate polynomials. 

We denote by (M[x] t )* the space of linear functionals on the set R[x] t of polynomials with degree 
at most t and use the notion of moment matrix M S (L) := (L(x Q x^)) (indexed by monomials of 
degree at most s) for L G (R[x]2 a )*- (See Section 2 for more definitions.) Say we want to compute 
the (finite) real variety Vk(I) of an ideal I given by a set of generators hi, . . . , h m G R[x] with 
maximum degree D. A common step in both methods is to compute a maximum rank moment 
matrix My t / 2 \ (L), where L G (M[x] f )* vanishes on the set Tit of all prolongations up to degree t of 
the polynomials hj ; this step is carried out with a numerical algorithm for scmidefinite optimization. 
From that point on both methods use distinct strategies. In the moment-matrix method one checks 
whether the rank condition: rankM s (L) = rankM s _i (L) holds for some D < s < [t/2\; if so, then 
one can conclude that \/l is generated by the polynomials in the kernel of M S (L) and extract 
Vr(I); if not, iterate with t + 1. In the prolongation-projection algorithm, one considers Gt, the set 
obtained by adding to Tit prolongations of the polynomials in the kernel of My t / 2 ](L), its border 
C? t + := G t Ui XiGt, as well as the set Gt °f linear functionals on IR[x] t vanishing on Gt, and its 
projections ir s (G^) on various degrees s < t. We give conditions on the dimension of these linear 
subspaces ensuring the computation of the real variety Vr(I) and generators for the real radical 
ideal \fl. Namely, if dim 7r s = dim it s -i(Gt~) = dini7r s ((C/ f + )- L ) holds for some D < s < t, 
then one can compute an ideal J nested between / and v^T so that Vjgi(I) = Vr(J), with equality 
J = \fl if dim7r s (C/ t - L ) = |Vr(/)|; if not, iterate with t+1. 

Both algorithms are tailored to finding real roots and terminate assuming Vr(7) finite (while 
Vc(I) could be infinite). However, the order t at which the dimension condition holds is at most 
the order at which the rank condition holds. Hence the prolongation-projection algorithm ter- 
minates earlier than the moment-matrix method, which often permits to save some scmidefinite 
optimization step with a larger moment matrix (as shown on a few examples in Section 6). 

Contents of the paper. Section 2 provides some basic background on polynomial ideals and 
moment matrices whereas Section 3 presents the basic principles behind the prolongation-projection 
method and Theorem 4, our main result, provides a new stopping criterion for the computation 
of Vr(J). Section 4 relates the prolongation-projection algorithm to the moment-matrix method 
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of [12]. In particular, Proposition 12 shows that the rank condition used as stopping criterion in 
the moment-matrix method is equivalent to a strong version of the new stopping criterion; as a 
consequence the new criterion is satisfied at least as early as the rank condition (Corollary 13). 
Section 5 contains a detailed description of the algorithm whose behavior is illustrated on a few 
examples in Section 6. 

2. PRELIMINARIES 

2.1. Polynomial ideals and varieties. We briefly introduce some notation and preliminaries 
for polynomials used throughout the paper and refer e.g. to [4] and [3] for more details. 

Throughout R[x] := R[a;i, . . . , x n ] is the ring of real polynomials in the n variables x = 
(x\, . . . , x n ) and R[x] t is the subspace of polynomials of degree at most t 6 N. For a G N", 
x" = x" 1 • • • x" n is the monomial with exponent a and degree |a| = J2i a i- For an integer t > 0, 
the set N™ = {a G N n \ \a\ < t} corresponds to the set of monomials of degree at most t, and 
T™ = {x Q | a G N"}, T" = {x Q | a G N™} denote the set of all monomials and of all monomials 
of degree at most t, respectively. Given S C R[x], set XiS :— {xip \ p G S}. The set 

S + := 5U.ti5U ...Ux n S 

denotes the one degree prolongation of S and. for B C T", dB := B + \ B is called the set of 
border monomials of B. The set B is said to be connected to 1 if any m G B can be written as 
m = mi ■ ■ ■ mfc with m\ = 1 and m\ ■ ■ ■ mn € B for all h = 1, . . . , k. For instance, B is connected 
to 1 if it is closed under taking divisions, i.e. racB and ml divides m implies m! G B. 

Given h\, . . . , h m G R[x], / = (h\, . . . , h m ) is the ideal generated by h\, . . . , h m , its algebraic 
variety is 

V C (I) := {v G C" | hj(v) = Vj = 1, . . . ,m)} 

and its real variety is Vr(I) := M™ n Vc(/). The ideal J is zero-dimensional when Vc(I) is finite. 
The vanishing ideal of a set V C C" is the ideal 

/(F) := {/ G R[x] | f(v) = Vv G V}. 

The Real Nullstellensatz (see e.g. [4, §4.1]) asserts that /(Vr(/)) coincides with ^7, the real radical 
of I defined as 

VI := \p G R[x] | p 2m + ^ 9 ? G / for some q. } G R[x] , m G N \ {0}}. 

3 

Given a vector space A on R, its dual vector space is the space A* = Hom(yl, R) consisting of 
all linear functionals from A to R. Given B C A, set B 1 - := {L e A* \ L(b) = V6 G B}, and 
Span R (B) := {£™ i Kh \ X t G R, b l G £?}. Then Span R (B) C (B^)^, with equality when A is 
finite dimensional. 

For an ideal I C R[x], the space V[I] := I 1 ^ ={Le (R[x])* | L(p) = Vp G /}, considered 
e.g. by Stetter [26], is isomorphic to (R[x]//)* and 2?[i] = / when I is zero-dimensional. Recall 
that I is zero-dimensional precisely when dimR[x]/J < oo, and |Vc(/)| < dimR[x]/7 with equality 
precisely when I = I(Vc(I))- 

The canonical basis of R[x] is the monomial set T™, with T> n := {d a |G N™} as corresponding 
dual basis for (R[x])*, where 

da(p) = TT^! { dx/ldx^ ) (0) fo ^ eR M- 

Thus any L G (R[x])* can be written in the form L = ^ Q y a d a (for some y G R N ). 
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By restricting its domain to R[x] s , any linear form L £ (R[x])* gives a linear form n s (L) in 
(R[x] s )*. Throughout we let ir s denote this projection from (R[x])* (or from (R[x] t )* for any 
t > s) onto (R[x] s )*. 

Given a zero-dimensional ideal 7 C R[x], a well known method for computing Vc(I) is the 
so-called eigenvalue method which relics on the following theorem relating the eigenvalues of the 
multiplication operators in R[x]/7 to the points in Vc(7). See e.g. [3, Chapter 2§4] . 

Theorem 1. Let I be a zero- dimensional ideal in R[x] and h £ R[x]. The eigenvalues of the 
multiplication operator 

m h : R[x]/7 — > R[x]/7 

p mod 7 i ► ph mod 7 

are the evaluations h(v) of the polynomial h at the points v £ Vc(7). Moreover, given a basis B 
o/R[x]/7, the eigenvectors of the matrix of the adjoint operator of with respect to B are the 
vectors (b(v)) beB £ IK. I ^ I (for all v £ V C {I)). 

The extraction of the roots via the eigenvalues of the multiplication operators requires the 
knowledge of a basis of R[x]/7 and an algorithm for reducing a polynomial p e K[x] modulo the 
ideal 7 in order to construct the multiplication matrices. Algorithms using Grobner bases can be 
used to perform this reduction by implementing a polynomial division algorithm (see [4, Chapter 
1]) or, as we will do in this paper, generalized normal form algorithms using border bases (see [13], 
[21], [26] for details). 

2.2. Moment matrices. Given L £ (R[x])*, let Ql denote the quadratic form on R[x] defined 
by Ql{p) '■= L{p 2 ) for p £ R[x]. Ql is said to be positive semidefinite, written as Ql X 0, if 
Ql(p) > for all p £ R[x]. Let M(L) denote the matrix associated with in the canonical 
monomial basis of R[x], with (a,/3)-entry 7(x c *x' 3 ) for a,f3 £ N n , so that 

Ql(p)= V a PpL{x. a ^) = vec(p) T Af (7,)vec(p), 

where vec(p) is the vector of coefficients of p in the monomial basis T". Then Ql X if and 
only if the matrix M{L) is positive semidefinite. For a polynomial p £ R[x], p £ KerQ^ (i.e. 
Qlip) = and so L(pq) = for all q £ R[x]) if and only if M(L)vec(p) = 0. Thus we may identify 
Ker M(L) with a subset of R[x], namely we say that a polynomial p £ R[x] lies in Ker M(L) if 
M(L)vec(p) = 0. Then KerA7(7) is an ideal in R[x], which is real radical when M(L) X (cf. 
[15], [17]). For an integer s > 0, M S {L) denotes the principal submatrix of M{L) indexed by N™. 
Then, in the canonical basis of R[x] s , M S (L) is the matrix of the restriction of Ql to R[x] s , and 
Ker M S (L) can be viewed as a subset of R[x] s . It follows from an elementary property of positive 
semidefinite matrices that 

(1) M t (L) X =*> Ker M t (L) n R[x] s = Ker M S (L) for 1 < s < t, 

(2) M t (L), Mt(L') y Ker M t (L + L') = Ker M t (L) n Ker M t (L'). 

We now recall some results about moment matrices which played a central role in our previous 
work [12] and are used here again. 

Theorem 2. [5] Let L £ (R[x] 2s )*. 7/rankM s (7) = rankM s _i(L), then there exists (a unique) 
L £ (R[x])* such that ir 2s {L) = L, rankM(L) = rankM s (7) ; and Ker M(L) = (KcrA7 s (7)). 

Theorem 3. (cf. [12, 15}) Let L £ (R[x])*. If M(L) h and rankA7(i) = rankM s _i(7), then 
Ker A7(7) = (Ker M S (L)) is a zero- dimensional real radical ideal and |Vc(Ker M (L))\ = ranki\7(7). 
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3. Basic principles for the prolongation-projection algorithm 

We present here the results underlying the prolongation-projection algorithm for computing 
Vk(I), K. = R, C. The basic techniques behind this section originally stem from the treatment 
of partial differential equations. Zharkov et al. [30, 31] were the first to apply these techniques 
to polynomial ideals. Section 3.1 contains the main result (Theorem 4). The complex case is 
inspired from [32] and was treated in [13]. The real case goes along the same lines, so we only 
give a brief sketch of proof in Section 3.2. In Section 3.3 we indicate a natural choice for the 
polynomial system Q involved in Theorem 4, which is based on ideas of [12] and will be used in 
the prolongation-projection algorithm. 

3.1. New stopping criterion based on prolongation/projection dimension conditions. 

We state the main result on which the prolongation-projection algorithm is based. We give a 
unified formulation for both complcx/rcal cases. 

Theorem 4. Let I = (hi, . . . , h m ) be an ideal in R[x], D = max., deg(hj) and s,t be integers with 
1 < s < t. Let Q C K[x]j satisfying hi, ... , h m € Q and SCI (resp., Q C \fl). If &\vci'K 3 (Q^) = 
then Vc(-0 = (resp., Vr(7) = 0). Assume now that s> D and 

(3a) dimTT s (g ± ) = dxmir s -i(Q ), 

(36) dim tt^) =dimn s ((g+) ± ). 

Then there exists a set B C T"_ x closed under taking divisions (and thus connected to 1) for which 
the following direct sum decomposition holds: 

(4) R[x] s = Span R (6) (R[x] s n Span R (£)). 

Let B C T"_j be any set connected to 1 for which (4) holds, let (p be the projection from R[x] s onto 
Span R (S) along R[x] s n Span K (C/), and let Fq := {ip(m) \ m G dB}, J := (Fq). Then B is a basis 
o/R[x]/J and Fq is a border basis of J . Moreover: 

• IfQQI then J = 1. 

• IfQQ\fl then 

Vr(I) = V C (J) n R" ; J n R[x] s = Span E (^) n R[x] s ; n s (V[J]) = tt,^), 
and in addition, J = \fl i/dim7r s (C/" L ) = |Vr(/)|. 
This result is proved in [13] in the case when Q = H. t C /, where 

(5) Ht := {^hj | \a\ + deg(^) < t, j = 1, . . . , m} 

consists of all prolongations to degree t of the generators hj of /. Note however that in [13] we did 
not prove the existence of B closed under taking divisions; we include a proof in Section 3.2 below. 

The proof for arbitrary Q C I is identical to the case Q — Tit- In the case Q C v 7 /, the proof 1 
is essentially analogous (except for the last claim J = \/l which is specific to the real case). We 
give a brief sketch of proof in the next section, since this enables us to point out the impact of the 
various assumptions and, moreover, some technical details are needed later in the presentation. 



Note that if we would apply the previous result to the ideal J := (I U Q) and the set Q, then we would reach 
the desired conclusion, but under the stronger assumption s > max(D, D'), where D' is the maximum degree of a 
generating set for Q. 
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3.2. Sketch of proof for Theorem 4. We begin with a lemma used to show the existence of B 
closed by division in Theorem 4. 

Lemma 5. Let Y be a matrix whose columns are indexed by T™. Assume 

(6) VAeR^l \>Y a = 0=> Yl X aY Xta =Q, 

where Y a denotes the a-th column ofY. Then there exists B C T" which is closed under taking 
divisions and indexes a maximum linearly independent set of columns ofY. 

Proof. Order the monomials in T™ according to a total degree monomial ordering ^. Let 
index a maximum linearly independent set of columns of Y , which is constructed using the greedy 
algorithm applied to the ordering ^ of the columns. Then, setting B m := {m' £ B \ m' -< m}, 
m £ B precisely when B m U {m} indexes a linearly independent set of columns of Y. We claim 
that B is closed under taking divisions. For this assume m £ B and m = i,mi with m\ ^ B. As 
mi $ B, we deduce that 

Y mi = X °Y* for 

some scalars A Q . 

For a £ B mi , a -< mi implies X{a -< Xiini = to, i.e., xia £ 2? m . Applying (6) we deduce that 

Y m ^ ^ ^aY Xi a: 

a€B mi 

which gives a linear dependency of Y m with the columns indexed by B m , contradicting m £ B. □ 

We now sketch the proof of Theorem 4. Set N := dim^-i^). If N = then Vfc(J) = 
(for otherwise the evaluation at v £ Vk{I) would give a nonzero element of it s -i{G^))- Let 
{Li, ...,Ljv} C for which {7T s _x(ii), . . . , 7T s _i(ijv)} is a basis of 7r s _i(^- L ). Let Y be the 
A x IT™^ matrix with (j,m)-th entry Lj{m) for j < N and to £ T^j. We verify that Y satisfies 
the condition (6) of Lemma 5 (replacing s by s— 1). For this note that X^aeT" ^a,Y a = if and only 
if p := J2 a eT"_ 2 £ (7r s _ 2 (^ ± )) ± = Span R (£) n R[x] s _ 2 and thus x t p £ Span R (£+) n R[x] s _i; 
in view of (36), this implies XiP £ Span R (C7) DM[x] s _i and thus X^aeT™ ^a,Y Xi a = 0. Thus we can 
apply Lemma 5: There exists a set B indexing a maximum linearly independent set of columns of 
Y which is closed by division. This amounts to having the direct sum decomposition: 

(7) R[x]._i = Span ffi (B) (Span R (£) fl R[x] s _ x )). 

As A = dimTT s (Q- L ), the set {ir s (Li), . . . ,it s (Ln)} is a basis of tt s (G ± ), and thus (4) holds. Set 
F := {to — (f(m) | to £ T™}. Obviously, Fq C F C Span R (C/ t ) n M[x] s . Moreover, one can verify (cf. 
[13]) that 

(8) Span R (F) = Span R (£) (1 R[x] s , 

(9) (F ) = (F), I C (F) ]£ 8 >D, 

(10) (p(xi(f(xjm)) = ip(xj(p(xim)) for m £ B and i = 1, . . . , n. 

Note that (36) is used to show (9)-(10). 

The ideal J := (i*b) satisfies 7 C J (by (9)) and J £ I or J £ \fl depending on the assumption 
on C?. As B is connected to 1 and we have the commutativity property (10), we can apply [18, 
Theorem 3.1] and deduce that B is a basis of M[x]/J. The inclusion: Span R (£) n R[x] s C J n 
R[x] s follows from (8)-(9), while the reverse inclusion follows from the fact that tp(p) = for 
all p £ J n R[x] s since B is a basis of R[x]/J. Thus Span R (C/) R[x] s = Jfl R[x] s , implying 
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tt s (Q ± ) = (Jnlfx],,) 1 . The inclusion w s (J ) C (JflM[x] s ) is obvious, and the reverse inclusion 
follows from (7T S (J )) C (J ) _L nR[x] s = Jnl[x] s , since J is zero-dimensional. Hence tt s (Q ± ) = 
n s (J ± ) = n s {V[J]). Finally note that 



Hence, if dim7r s (CJ- L ) = |Vr(/)|, then equality holds throughout, which implies that J is real radical 
and thus J = \fl. This concludes the proof of Theorem 4. 

Remark 6. We indicate here what happens if we weaken some assumptions in Theorem 4- 

(i) The condition s > D is used only in (9) to show I C (F). Hence if we omit the condition s > D 
in Theorem 4, then we get the same conclusion except that we cannot claim I C J. 

(ii) Consider now the case where we assume only that (3a) holds (and not (3b)). As we use (3b) to 
show the existence of B connected to 1 and to prove (9)- (10), we cannot prove the commutativity 
property (10), neither the equality (F) = (Fq). Nevertheless, what we can do is test whether B is 
connected to 1 and whether (10) holds. If this is the case, then we can conclude that B is a basis 
o/R[x]/J where J = (Fq) C ^7 and extract the variety Vc(J) which satisfies Vr(I) C Vc(J) (11" 
and |Vc(J)| < dimM[x]/J = \B\. Then it suffices to sort Vr(I) out of Vc(J). The additional 
information that condition (3b) gives us is the guarantee that the commutativity property (10) 
holds and the equality J — (F), thus implying J D I and Vr(I) = Vc(J) PI W l if s > D. 

3.3. A concrete choice for the polynomial system Q in Theorem 4. For the task of com- 
puting Vc(I), one can choose as indicated in [13] the set Q = Tit from (5) and thus consider the 
linear subspace JCt := Ti.^ of (]R[x] t )*. For the task of computing Vr(I), as inspired by [12], we 
augment Ti t with a set Wt of polynomials in \/l obtained from the kernel of a suitable positive 
element in H^~. For this, consider the convex cone 



consisting of the elements of K-t that are positive, i.e. satisfy L(p 2 ) > whenever deg(p 2 ) < t. 
Generic elements of ICt^y (defined in Lemma 7 below) play a central role; geometrically these are 
the elements lying in the relative interior of the cone K-t^. 

Lemma 7. The following assertions are equivalent for L* G K-t.y- 

(i) rankM Lt/2 j (L*) = max LEKt > rankM Lt/2 j (L). 

(ii) rankM s (L*) = max Le/Ct , t rankM s (i) for all 1 < s < [t/2\ . 

(iii) Ker M S (L*) C Ker M S (L) for all L G K tj y and 1 < s < [t/2\. 

Then L* is said to be generic. 

Proof. Direct verification using (l)-(2). □ 

Hence any two generic elements L\,L2 G K-t,y have the same kernel, denoted by Aft (= 
KerM[ t / 2 j(Li) = KerM^/ 2 j (^2)), which satisfies 



dimTT.^) = \B\ = dimK[x]/J > |y c (J)| > \V R {I)\. 



Kt, h :={L&Ht\M vm {L)±Q}, 
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Therefore, W t C \fl. For the task of computing Vr(J), our choice for the set Q in Theorem 4 is 

(15) &:=WtUW t . 
Note also that 

(16) iCt, h c^nw^fWtUWt) 1 . 

In fact, as we now show, both sets in (16) have the same dimension, i.e. (hit U Wt)" 1 " is the smallest 
linear space containing the cone K-t^- 

Lemma 8. dim/C t > = <tim(H t U W t )^- 

Proof. Pick L* lying in the relative interior of JCt,y, i-c. L* is generic, and define 

T t := {L e (R[x] t )* | L* ± eL e for some e > 0}, 

the linear space consisting of all possible perturbations at L*. Then, dim/C t > = dim"P t . One 
can verify that there exists an e > such that L* ± eL £ fct,y if and only if L £ and 
KerM Lt /2j (L*) C Ker M|t/ 2 j (L) (cf. e.g. [8, Thm. 31.5.3]). As the latter condition is equivalent 
to L £ by (14), we find P t = (H t U Wt)" 1 , which concludes the proof. □ 

We conclude with a characterization of and of its dual space 2? [-^7], using the sets Qt from 
(15). 

Proposition 9. With Q t = H t U W t , ^1 = \J t Span R (£ t ) and V[yfl\ = f] t . 

Proof. The inclusion (J t Span R (C* t ) C \fl follows from (12). Next, for some order (t,s) we have 
\/l = (Ker M S (L*)). The proof, which relies on the existence of a finite basis for the ideal yl 
can be found in [12]. This fact, combined with Ker M S (L*) C Mt Q Span K (^t), implies the reverse 
inclusion C |J t Span R (£ t ). Now the equality y/l = \J t Span K (? t implies in turn 2?[v^T] = P| t Q^. 

□ 

When | Vr (I) \ < oo, the dual of the real radical ideal coincides in fact with the vector space 
spanned by the evaluations at all v £ Vr(J). Proposition 9 shows how to obtain it directly from 
the quadratic forms Ql (or its matrix representation My t / 2 \ (L)) for a generic L £ JCt,y without a 
priori knowledge of Vr(I). 

4. Links with the moment-matrix method 

In this section we explore the links with the moment-matrix method of [12] for finding Vr(7) as 
well as the real radical ideal v^T. We recall the main result of [12], underlying this method. 

Theorem 10. [12] Let L* be a generic element of K-t.y. Assume that 

(17) rankM s (L*) = rankM s _i(L*) 

for some D < s < [t/2\. Then (Ker M B (L*)) = \fl and any set B C T™_ 1 indexing a maximum 
linearly independent set of columns of M s _i(L*) is a basis o/l^xj/v 7 /. 
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4.1. Relating the rank condition and the prolongation-projection dimension condi- 
tions. We now present some links between the rank condition (17) and the conditions (3a)-(36). 
First we show that the condition (3 a) suffices to ensure that the rank condition (17) holds at some 
later order. 

PROPOSITION 11. Let l<s<t. If (3a) holds withQ := H t UW t , then va,nkM s (L) = rankM s _i(i) 
for all L G JC t+2 s,y- 

Proof. Let L G JCt+2s>- We show that rankM s (L) = rankM s _i(Z). For this, pick m,m' G 
T™. As in the proof of Theorem 4, (4) holds and thus we can write m = J2beB^ b ^ + /> wnere 
\b G R, / G Span R (C/), and B C T"_ x . (Note that (36) was not used to derive this.) Then, 
mm! = X)be8 A&m'6 + m' f. It suffices now to show that L(m! /) = 0. Indeed this will imply 
M(£) m / >m = L(mm') = ^2 b£ g\L(m'b) = X^beB A&M(L) m / i b, that is, the mth column of M(L) is 
a linear combination of its columns indexed by 6 G 23, thus giving the desired result. 

We now show that L(m'g) = for all g G H t U>Vt. By assumption, L G /C t+ 2s> C H£+ 2s nW t -+ 2s 
(recall (16)). If g € Ht, then m'g G W{+ s C Ht+2s and thus L(m'g) = 0. If <? G W*, then g = x"/i, 
where h G M and |a| < [i/ 2 J • Hence, m'g = m'x Q /i, where deg(m'x Q ) < s + [t/2j < L(2s + 1)/2\ 
and h € N't Q ■M+2s (by (11)), implying m'g G Wt+2s and thus L(m'g) = 0. □ 

We now show that the rank condition (17) is in fact equivalent to the following stronger version 
of the conditions (3a)-(36) with Q = Q t = Ti t U Wt- 

(18a) dim7T2 S (<?t ) = dim7r s -i(</ t ), 

(186) dim 7 r 2s (S i ± ) = dini7r 2s ((2+) ± ). 

PROPOSITION 12. Let L* be a generic element of JC t ,y and 1 < s < L*/2J . 

(i) Assume (17) holds. Then (18a) holds, and (18b) holds as well if s> D. 

(ii) Assume (18a)-(18b) hold. Then, (17) holds, the ideal J obtained in Theorem 4 is real 
radical and satisfies J = (Ker M S (L*)) C I(Vr(I)) and, given B C T"_i, B satisfies (7) 
if and only if B indexes a column basis of M S —±(L*). Furthermore, J=v r fifs>D. 

The proof being a bit technical is postponed to Section 4.2. An immediate consequence of Propo- 
sition 12 is that the rank condition at order (t, s) implies the prolongation-projection dimension 
conditions (3a)-(36) at the same order (i, s). 

Corollary 13. Assume D < s < |_t/2j and let g = g t =H t UW t . Then, 

(17) (18a)-(18b) =S> (3a)-(3b). 

Proof. Indeed, ir s (Q^) = 'n's{{Gt) ± follows directly from TT2 S (Gt~) = 7r 2s ((£7^ ■ ^ 

It is shown in [12] that the rank condition (17) holds at order (s, t) large enough with D < s < 
[_i/2_|. Hence the same holds for the conditions (18a)-(186) (and thus for (3a)-(36)), which will 
imply the termination of the prolongation-projection algorithm based on Theorem 4. 

4.2. Proof of Proposition 12. First we note that the rank condition (17) is in fact a property 
of the whole cone ICt.y and its superset = TC^~ fl W^. 

Lemma 14. If (17) holds for some generic L* G JC t .y, then (17) holds for all L* G . 
Proof. Let L G G± ■ Wc have 

(19) Ker M S (L*) = Ker M Lt/2J (L*) n R[x] s = Af t R R[x] s C Ker M Lt/2J (L) n R[x] s C Ker M S (L), 
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where the first equality holds by (1), the first inclusion holds by (14), and the second one holds 
since M S (L) is a principal submatrix of M^ t / 2 \(L). This implies directly that rankM s (L) = 
rankM s _i(Z). □ 

We now give the proof for Proposition 12. Let L* be a generic element of ICt.y. 

(i) Assume that (17) holds. First we show (18a), i.e. we show that dim7r 2s (^ ( L ) = dim7r s _ 1 (^ t L ). 
For this, consider the linear mapping 

As i/j is onto, it suffices to show that ip is one-to-one. For this assume 7r s _i(L) = for some 
L G Gt- We show that tt 2s (L) = 0, i.e. i(x 7 ) = for all I7I < 2s by induction on | -y | < 2s. The 
case I7I < s — 1 holds by assumption. Let s < | — y | < 2s and write 7 as 7 = a + f3 where \a\ = s 
and \j3\ < s. By Lemma 14, rankM s (L) = rankM s _i(L). Hence the ath column of M S (L) can be 
written as a linear combination of the columns indexed by T™^. This gives 

M s {L) p . a = £|5|< s _i X s M s (L) 0! s for some X s G R. As \/3 + 5\ < \-y\ - 1, we have M s (L) 0iS = 
L(yfi +5 ) — by the induction assumption, implying I/(x 7 ) = M s (L)p^ a = 0. 

We now assume moreover s > D. We show the inclusion ~K2s(Gt) Q 7T 2s((G^)' L ), which implies 
(186). Let L G Gt ■ As rankM s (L) = rankM s _i(L), we can apply Theorem 2 and deduce the 
existence of L G (R[x])* for which 7r 2s (L) = 7r 2s (L) and Ker A/(Z) = (Ker M S (L)). It suffices now to 
show that L G (G^) ± - We show a stronger result, namely that L G I(Vj&(I)) . As s > D, we know 
from Theorem 10 that I(V R (I)) = (KeiM s (L*)). Pick p G I(V m (I)) and write it as p = Y.i u Wh 
where ui G R[x] and gi G Ker M S (L*); we show that L(p) = 0. By (19), g\ G Kev M S (L) and thus, 
as M S (L) = M S (L), g t G KerM s (L). Therefore, p lies in (KerM s (i)) = KeriW(L), which gives 

Hp) = 0. 

(ii) Assume now that (18a)-(186) hold. Then, (3 a) -(3b) holds for the pair (t, 2s) (and G = Gt)- 
Although we do not assume 2s > D, the conclusion of Theorem 4 partially holds, as observed in 
Remark 6 (i). Namely, we can find an ideal J satisfying J C /(Vr(/)), J n R[x] 2s = Span R (C/ t ) n 
R[x] 2s , iT2s(D[J}) = iT2s(Gt~), and / C J if 2s > D. Moreover, there exists a set B C T"_ x which 
is a basis of R[x]/J and satisfies the following analogue of (4): 

(20) R[x] 2s = Span K (B) 8 (S P an R (£ t ) n R[x] 2s ). 

We show that rankM s (Z*) = rankM s _i(L*). As L* G G^~, there exists L G D[J] for which 
K2s(L*) = tt 2s (L). Thus M S (L*) = M S (L), and J C KerM(Z) since L G D[J]. It suffices to show 
that r&nkM s (L) = rankM s _i(L). For this, as in the proof of Proposition 11, pick m,m' G T". 
Using (20), we can write m = X^beB ^bb+f, where A& G R, / G Span R (£ t )nR[x] 2s C J C Ker M(L), 
so that L(m'm) = X^ee XbL(m'b)^ which gives the desired result: rankA/ s (L) = rankM s _i(L). 
Let Si,i3 2 C T™_ x , where B\ satisfies (7) and B2 indexes a column basis of M s _x(L*). Then 

(21) \Bx\ = dimR[x]/J < rankM s _i(i*) (= \B 2 \) 

since the columns of M S -\(L*) indexed by B\ are linearly independent (direct verification, using 
(7) and the fact that KerM s _i(L*) C Ker M [t/2] (L*) =Af t C Span R (£ t )). Moreover, 

(22) |B 2 | = rankM s _!(L*) < dim7r a _i(0 t x ) (= |Bi|). 

Indeed, as Span R (C/ t ) n R[x] s _i C J n R[x] s _i C KerM s _i(L) = KerM s _i(L*), we obtain 
Span R (£ t ) n Span R (S 2 ) = {0}, which implies \B 2 \ < dim(Span R (^ t ) n R^^^J^ = diiinr s ^i(G^). 
Hence, equality holds in (21) and (22). Therefore, B\ indexes a column basis of M s _i(L*), S 2 
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satisfies (7), and 

rankM s _i(L*) = dim7r a _i(#-) = dimR[x]/J. 
As J C Ker M(L), we deduce 

dim R[x] /Ker M(L) < dimR[x]/J. 

On the other hand, 

dimR[x]/J = rankM s _i(£*) = rankA/ s _i(L) < rankM(L) = dimR[x]/KerM(L). 

Hence equality holds throughout. In particular, J = Ker M{L) and rankAf(L) = rankM s _!(i). 
As M a -i(L) = M a -i(L*) h 0, we deduce that M{L) > and J = KerAf(L) = (KerM s (Z)) = 
(Ker M S {L*)) is a real radical ideal (using Theorem 3). Finally, if s > D, then J = (KerM s (L*)) = 
\fl by Theorem 10. This concludes the proof of Proposition 12. 

4.3. Two illustrative examples. We discuss two simple examples to illustrate the various no- 
tions just introduced and the role of moment matrices; the second one has infinitely many complex 
roots. 

Example 15. Let I = {x\,x\,x\x 2 ) C IR[xi,X2], considered in [13] as an example with a non- 
Gorenstein algebra R[x]/7. Any L £ ICt (t > 2 ) satisfies L(x. a ) = if \a\ > 2 and thus 



b 
c 




for some scalar s a,b,c, 



) 

im7T2(/C 2 ) 



dim7Ti(/C2) = dim7T2(/C3) = 3 and 



where entries are indexed by 1, x\, x%, . . . Hence, di 
the rank stabilizes at order (t, s) = (4, 2), i.e. rank/I^ (L*) = rankMi(L*) = 2 for generic L* £ /C4. 
When L € Kt,y, the condition M^ t / 2 \(L) y implies b = c = 0. Hence, for generic L* G /C2>? 
Ml '■= KerAfi(L*) is spanned by the polynomials X\ and x 2 , and the rank condition (17) holds 
at order (t,s) = (2,1), i.e. rankMi(L*) = rankMo(L*) = 1. As Span R (C/2) is spanned by the 
polynomials x\,x 2 ,x\,x\X2,x\, the conditions (18a)-(186) hold at the same order (t,s) = (2,1), 
i.e. dim7T2(^ 2 1 ) = dim7ro(^ 2 L ) = dim7T2((^ 2 f )" L ) = 1; as predicted by Proposition 12. 

Example 16. Consider the ideal I = {x\ + x\) C M.[x\,X2\ with Vr(I) = {0} and \Vc{I)\ = 00. 
As dim7r s (/Ct) = dini7r s _i(/Ct) + 2 for any t > s > 2, the conditions (Sa)-(Sb) never hold in the 
case Q = Jit- On the other hand, any L € IC2,y satisfies L(x\) = L{x\) = 0, which follows from 
L(x\ + X2) = combined with M\{L) > 0, giving L{x\), L(x 2 ) > 0. Moreover, L(x\) = L(x 2 ) = 
L(x\X2) = 0. Thus A/2 is spanned by the polynomials X\ and x 2 , and the conditions (17) and 
(18a)-(18b) hold at order (t,s) = (2, 1). 

Examples 18 and 20 in Section 6 are cases where the prolongation-projection method terminates 
earlier than the moment-matrix method. 



5. A PROLONGATION-PROJECTION ALGORITHM 

Let us now give a brief description of our algorithm for computing Vk(I) (K = R, C) based on 
the results of the previous section. A simple adjustment in the proposed prolongation-projection 
algorithm allows the computation of all complex vs. real roots. The general structure is shown in 
Algorithm 1. If I is an ideal given by a set of generators and |Vk(-0| < 00, this algorithm computes 
the multiplication matrices in R[x]/J, which thus allows the immediate computation of Vk(J) (by 
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Theorem 1), where J is a zero-dimensional ideal satisfying J = / if K = C and I C J C \fl if 
K = R, so that Vk(J) = Vk(I). We then comment on the key steps involved in the algorithm. 



Algorithm 1 Unified prolongation-projection algorithm for computing Vk{I): 
Require: A set {hi, ... , h m } of generators of I and t > D. 

Ensure: The multiplication matrices in R[x]/J, where J = I if K = C and I C J C \fl if K = R, 
thus enabling the computation of Vk(I). 
1: Compute the matrix representation Gt of Gt and of Gt ■ 
2: Compute KerG t and KerG+. 

3: Compute dim 7r s (KerGt) (= diim: s ((Gt)~ L )) and dim n s (Ker Gf ) (= dim7r s ((tJ t + )- L )) for s < t. 
4: Check if (3o)-(36) holds for some D<s< [t/2\. 
5: if yes then 

6: return a basis B C R[x] s _i connected to 1 and satisfying (7), and the multiplication 

matrices Xi in R[x]/J represented in the basis B. 
7: else 

8: Iterate (go to 1) replacing t by t + 1. 
9: end if 

Remark 17. ifere, 5* = 7it (see (5)) for the task of computing Vc(I), and Gt = Ht U W* fsee 
f or th- e t as k °f computing Vr(/). S'ee 6e/ow /or details about the matrix representations Gt 
and G^ . 



Characterizing Gt and Gt v ^ a ^ e niatrix Gt- In the real case, the set Gt is defined as Gt = TitUWf 
where Wt is the linear space defined in (13). As we are interested in the orthogonal space Gt, it 
suffices to compute a basis Ct of the linear space Aft and to define the set 

(23) St:={* a g\\a\<\t/2\,geCt}. 

Then, Aft = Span R (C t ), Wt = Span R (6>t), and Gt~ = (Ht U $t) ± - Let St (resp., H t ) be the matrix 
with columns indexed by T" and whose rows are the coefficient vectors of the polynomials in St 
(resp., in Tit). In the case K = C, the set Gt = 'Ht is represented by the matrix Gt ■= H t and, in 
the case K = M, the set Gt = ~Ht U Wt is represented by the matrix 



Gt 



Ht 
St 



Then the vectors in KerGt are precisely the coordinate vectors in the canonical basis of (R[x] t )* 
of the linear forms in Q^~, i.e. 

(24) Le^^(L(x Q )) H <tGKerGt. 

Analogously, Gf is the matrix representation of (HtU St) + , so that {Gt) 1 ' corresponds to KerG^ . 

To compute the space Aft we need a generic element L* € fct,y- How to find such a generic 
element has been discussed in detail in [12, Section 4.4.1]. Let us only mention here that this task 
can be performed numerically using a standard semidefinite programming solver implementing a 
self-dual embedding strategy, see e.g. [7, Chapter 4]. For our computations we use the SDP solver 
SeDuMi [27]. 
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Computing tt s (Q^~) and its dimension. As shown in (24), the dual space Gt~ can he characterized in 
the canonical dual basis as the kernel of the matrix Gt, see e.g. [32] for details using an algorithm 
based on singular value decomposition. Faster implementations can be obtained e.g. using Gauss 
elimination. Once we have a basis of KcrGf, denoted say by {z\, . . . , zm}, then, for any s < t, we 
construct the matrix Z s whose rows are the vectors ix s (zx), . . . , tt s (zm), the projections onto R" of 
Z\, ... , zm- Then dimir s (Gi~) = dini7r s (KerGt) is equal to the rank of the matrix Z s . 

Extracting solutions. In order to extract the variety Vk(-0, we apply Theorem 1 which thus requires 
a basis B of the quotient space and the corresponding multiplication matrices. In the setting of 
Theorem 4, rankZ s = rankZ s _i —: N and B is chosen such that B C T"_ x indexes N linearly 
independent columns of Z s -i. A first possibility to construct B is to use a greedy algorithm 
as explained in the proof of Lemma 5. Another possibility is to use Gauss- Jordan elimination 
with partial pivoting on Z s (see [10]) such that each column corresponding to a monomial of 
degree s is expressed as a linear combination of N monomials of degree at most s — 1. The 
pivot variables form a set B C T"_i indexing a maximum set of linearly independent columns 
of Z s and their corresponding monomials serve as a (monomial) basis B of the quotient space 
(provided B is connected to 1). The reduced row echelon form of Z s , interpreted as coefficient 
vector for some polynomials, gives the desired rewriting family, which thus enables the construction 
of multiplication matrices and provides a border (or Grobner) basis (cf. [12] for details). 

A second alternative proposed in [32] is to use singular value decomposition once more to obtain 
a basis of Ker Z s and therefore a polynomial basis B for the quotient ring (see [32] for details). All 
examples presented in the next section are computed using singular value decomposition. 



We now illustrate the prolongation-projection algorithm on some simple examples. The algo- 
rithm has been implemented in Matlab using the Yalmip toolbox [16]. For the real-root prolongation- 
projection algorithm, we show the dimensions of ir s (G^) and Ks((Gt)' L )> the projections of the or- 
thogonal complement of the set Qt = Ht U Wt and of its one degree prolongation. For comparison, 
we also sometimes show the dimension table for the complex-root version of this algorithm, and we 
show the values rankM s (i*) (s < |_£ /2j ) for a generic element L* € fct,y used in the real moment- 
matrix method. To illustrate the potential savings, and at the same time facilitate the comparison 
between the various methods, we sometimes give more data than needed for the real root com- 
putation (then displayed in color gray). We also provide the extracted roots v £ Vk(-0 and, as a 
measure of accuracy, the maximum evaluation e(v) = maxj |/ij(w)| taken over all input polynomials 
hj at the extracted root v, as well as the commutativity error c(X) := max" 3=1 abs(XiXj — XjXi) 
of the computed multiplication matrices Xi. 

Example 18. Consider the ideal I = (hi, /i2 ) h^) C R[xi, X2, X3], where 



with D = 3 ; |Vc(/)| = 8 and |Vr(/)| = 2 ; taken from [3, Ex. 4, p. 57]. We illustrate and compare 
the various algorithms on this example. 

Table 1 shows the dimensions of the sets 7r s (/C t ) for various prolongation-projection orders (t, s). 
Note that the conditions (3a)-(3b) hold at order (t,s) = (6,3), i.e. 



6. Numerical Examples 



in 

h 3 



x\ - 2xixz + 5 , 

Xxxl + X2X3 + 1 , 
3^2 — 8X1X3 , 



dim TT3(ICe) — diimr2{ICe) = dim ^3(^7). 
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s = 





1 


2 


3 


A 


5 


6 


7 


8 


9 


dim7r s (/C 3 ) 


1 


4 


8 


11 














dim.7r s (/C4) 




4 


8 


10 


12 












dim7r s (/C 5 ) 


1 


4 


8 


9 


10 


12 










dim7r s (/C 6 ) 


1 


4 


8 


8 














dim7r s (/C7) 


1 


4 


8 


8 














dmnr s (JCs) 






















dim7r s (/Cg) 




4 



















Table 1. Dimension table for w s (]Ct) in Example 18. 



With the complex-root prolongation-projection algorithm we can compute the following eight com- 
plex roots: 



Vi 

t'2 
V3 
''4 
V5 
V6 
V 7 
V$ 



-1.10 -2.88 -2.82 ] , 

0.0767 + 2. 243i 0.461 + 0.497i 0.0764 + 0.00834i ] , 
0.0767- 2. 243i 0.461 - 0.497i 0.0764 - 0.00834i ] , 
-0.0815 - 0.931i 2.35 + 0.0431i -0.274 + 2.209i ] , 
-0.0815 + 0.931i 2.35- 0.0431* -0.274 - 2.20z ] , 
0.0725 + 2. 24i -0.466 - 0.464i 0.0724 + 0.00210i ] , 
0.0725- 2. 24i -0.466 + 0.464i 0.0724 - 0.00210i ] , 
0.966 -2.81 3.07 ] , 



with a maximum error of max^ e{v{) < 8e-13 and commutativity error c(X) < 6e-13. 



s = 





1 


2 


3 


4 


5 


6 


7 


dini7r s (^) 


1 


4 


8 


11 










dim 7 r s ((e 3 + ) ± ) 


1 


4 


8 


10 


12 








dim -K s {Gl) 


1 


4 


8 


10 


12 








dim^((5+)- L ) 


1 


4 


8 


9 


10 


12 






dim n s (Gs) 


1 


2 


2 


2 


3 


5 








1 


2 


2 


2 


3 


4 


6 




dim7r s (^6") 


1 




2 


2 


2 


2 


3 




dimir s ((G+) ± ) 


1 




2 


2 


2 


2 


2 




Dimension table for 






\ and 


*s{{Gt) 


x ) in 


E: 



Table 2 shows the dimensions of the sets tt s (G^) and ^ s ((Gt')' i ~) f or various prolongation- 
projection orders (t,s). Note that the conditions (3a)-(36) hold at order (t,s) = (5,2), i.e. 

dim7r 2 (C/g L ) = dim iti{G^) = dim 7r 2 ((C? 5 h ) ± ). 

With the real-root prolongation-projection algorithm we can extract the two real solutions: 

V! = [ -1.101 -2.878 -2.821 ] , 

v 2 = [ 0.966 -2.813 3.072 ] , 
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with maxi e(v{) < 2e-8 and commutativity error c[X) < 3.3e-9. Note that, since 2 = s < D = 3, we 
cannot directly apply Theorem 4 to claim Vr(7) = Vc{ J)nR™. Instead, as indicated in Remark 6 (i), 
we can only claim Vc(J) H ffi™ 2 Vr(7). However, equality can be verified by evaluating the input 
polynomials hj at the points v G Vc(J) flffi™. Anyway, one can also observe that the conditions 
(3a)-(36) hold at order (t,s) = (5,3), in which case one can directly conclude Vr(I) = Vc(J)nR n . 
Finally, we can even conclude J = ^fl since dim n s (Q^~ ) = |Vr(/)| (using the last claim in Theorem 

4)- 

The ranks of the moment matrices involved in the computation are shown in Table 3. Observe 
that the rank condition (17) holds at order (t,s) — (6,2), i.e. 

rankA/2(£*) = rankMi(L*) for generic L* G /C$ t y. 

(To be precise, as 2 = s < D = 3, we use [12, Prop. 4.1] and check whether the extracted roots 
belong to Vr(7) afterwards.) 





s = 


s = l s = 2 s = 3 


t = 


3 


1 


4 


t = 


4 


1 


4 8 


t = 


■5 


1 


2 8 


t = 


6 


1 


2 2 



Table 3. Showing rankM s (L*) for generic L* G JC t .y in Example 18. 



In this small example, we see that we can improve efficiency over the general complex-root 
algorithm if we are only interested in computing the real roots. Indeed the prolongation-projection 
algorithm terminates at order (t, s) = (5, 2) in the real case while it terminates at order (6, 3) in 
the complex case, however at the price of solving an SDP in the real case. Moreover, compared to 
the real-root moment-matrix algorithm of [12], we save the computation of the last moment matrix 
M 3 (L*) forL* G /C 6 >. 

Modifying the above example by replacing each polynomial hi by hi ■ (1 + x\) yields an example 
with a positive dimensional complex variety, while the real variety is unchanged. The proposed 
algorithm still converges, this time at order (t, s) = (7, 2) and allows the extraction of the two real 
roots. 

Example 19. Consider the ideal I = (hi, /12, /13) C M.[x±, X2], where 

h\ = x^x\ -\- 3x^ — X'2 — 

h% = x\x2 — 2xf, 

I13 = 2x\x\ — x\ — 2x\ + x\, 

and D — 5, taken from [3, p. 40]. The corresponding variety consists of two (real) points, one of 
which has multiplicity 8. 

Table 4 shows the dimensions of the projections of the sets Gt~ an d [Qt)^ ■ The conditions (3a) -(3&) 
hold at order (t, s) = (6, s) with 2 < s < 5, i.e. 

dim7r s (C/^) = dim7T s _i(</ 6 ) = dim7T s ((CJ 6 h ) J -) for 2 < s < 5, 

the conditions (18a)-(186) hold at order (t,s) = (6,2), i.e. 

dirndl (£7^-) = dim7r 4 (CJ^) = dim7r 4 ((£J 6 h ) ± ), 
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s = 





1 


2 


3 


4 


5 


6 


7 




1 


3 


5 


6 


8 


10 






dim7T s ((^) x ) 


1 


3 


5 


6 


6 


8 


10 




dim7r s (^) 


1 


2 


2 


2 


2 


2 


4 




dim7r s ((e+) x ) 


1 


2 


2 


2 


2 


2 


2 


4 



Table 4. Dimension table for n s (Gt~) and n s ((Gt~)~ L ) in Example 19. 
and the extracted roots are 

vi = [ -Q.VTe-6 1.10e-5 ] 
v 2 = [ 0.9988 1.9998 ] 

with an accuracy of e(v\) < 2e-10 and e(t>2) < 4e- 3 and maximum commutativity error c{X) < 3e-5. 
The ranks of the moment matrices involved in the computations are shown in Table 5. As predicted 
by Proposition 12, condition (17) holds at order (t,s) = (6,2), i.e. 

ranki\/2(X*) = rankMi(L*) for generic L* £ K,§ t y. 

Moreover, the returned ideal J satisfies J = (KerMi(L*)) = Table 6 shows the dimensions of 





s = 


s = 1 


s=2 s=3 


t = 5 


1 


3 


5 


t = 6 


1 


2 


2 4 



Table 5. Showing rankM s (L*) for generic L* £ K ti y in Example 19. 

the projections n s (JCt) for the complex-root prolongation-projection algorithm. The conditions (3a)- 
(3b) are satisfied at order (t,s) = (7,5), allowing (in principle) to extract the two roots with 
their corresponding multiplicities. The appearance of multiple roots requires a careful choice of the 
extraction procedure using multiplication operators. We employ the approach described in [6] using 
reordered Schur factorization. At order (t, s) — (7,5) , numerical problems prevent a successful 
extraction despite this algorithm. However, at order (t, s) = (8,5), the multiplication matrices (on 
which the reordered Schur factorization method is applied) have a commutativity error of c(X) < 
6.25e-16. Thus, we can extract the root 

v=[l 2] 

with accuracy e(v) < 1.38e-14 and the 8-fold root at the origin with an even higher accuracy of 
e(vi) < 1.75e-32. 

Note that the real version of this algorithm, working directly with the real radical of the ideal, 
does not require these considerations as it eliminates multiplicities. 

Example 20. This example is taken from [29] and represents a Gaussian quadrature formula with 
two weights and two knots, namely, I = (hi, . . . , h^), where 

h% = xi + X2 — 2 , 
h 2 = X1X3 + X2X4 , 

h 3 = xix\ + x 2 x\ - I , 
hi = xix 3 + X2X4, 
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s = 





1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


dim7r s (/C 5 ) 


1 


3 


6 


8 


11 


13 












dini7r s (/C 6 ) 


1 


3 


6 


8 


9 


11 


13 










dim 7r s (/C7) 


1 


3 


6 


8 


9 


9 












dmi7r s (/C 8 ) 


1 


3 


6 


8 


9 


9 












dim7r s (/Cg) 
























diiii7r s (/Ci ) 

























Table 6. Dimension tabic for ir s {Kt) in Example 19. 



with D = 4 and |Vr(/)| = |Vc(/)| = 2. Table 7 shows the dimensions for the projections of the sets 
and (£7 t + )~ L and Table 8 shows the ranks of the moment matrices M S (L*) for generic L* G /Ct>. 
The conditions (3a) -(36) hold at order (t,s) = (5,2) and the extracted roots are 

v x = [ 1 1 -0.5774 0.5774 ] 
v 2 =[l 1 0.5774 -0.5774 ] . 

with an accuracy of e(i>i) < 2e-ll and £(^2) < 2e-ll and maximum commutativity error c{X) < 
4e-14- Here again the algorithm returns the ideal J = since dim7r 2 (t/g L ) = |Vr(J)| = 2. On 
the other hand, the moment-matrix algorithm of [12] terminates at order (t, s) = (6,2), thus later 
than the prolongation-projection algorithm. 



s = 





1 


2 


3 


4 


5 


6 


7 


dim7r s (CJ^) 


1 


3 


7 


11 


20 








dim7r s ((5+) ± ) 


1 


3 


4 


8 


12 


23 








1 


2 


2 


2 


5 


16 






dim7r s ((£+H 


1 


2 


2 


2 


5 


9 


22 




dini7r s (^) 


1 


2 


2 


2 


2 


16 


18 




dim7r s ((£+H 


1 





2 


2 


2 


2 


2 


2 



Table 7. Dimension table for tt s (Q^) and 7r s ((C? f f )- L ) in Example 20. 





s = 


s = 1 


s=2 s = 3 


t 


= 4 


1 


4 


9 


t 


= 5 


1 


2 


5 


t 


= 6 









Table 8. Showing rankM s (L*) for generic L* G K. t ,y in Example 20. 



Example 21. The following 6-dimensional system is taken from 
http : //www .mat .univie . ac . at/~neum/glopt/ coconut /Benchmark/Library3/katsura5 .mod 
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and is known under the name Katsura 5: 

hi = 2x\ + 2x\ + 2x l + 2x i + 2x i +x\-xi, 
h 2 = x 6 x 5 + X5X4 + 2x4X3 + 2x 3 x 2 + 2x 2 x\ — x 2 , 
h 3 = 2x 6 X4 + 2x5X3 + 2x 4 x 2 + x\ + 2x^xi — xz , 
hi = 2x^X3 + 2x5X2 + 2x3X2 + 2x4X1 — X4 , 

= x\ + 2xqX\ + 2x5X1 + 2x4X1 — X5 , 
he = 2xq + 2x5 + 2^4 + 2x3 + 2x2 + x\ — 1, 

with D = 2, |Vc(-0| = 32, and |Vr(/)| = 12. The projection dimensions are shown in Table 9. 



s = 





1 


2 


3 


4 


5 


6 


7 


dim tt s (Q 2 ) 


1 


6 


16 












dim7T s ((S 2 + ) x ) 


1 


6 


16 


26 










dim7r s (^) 


1 


6 


16 


26 










dim7r s ((^) x ) 


1 


6 


16 


26 


31 








dini7r s (^) 


1 


6 


16 


26 


31 








dim^((^ 4 +)- L ) 


1 


6 


16 


26 


31 


32 






dim Ti s {Gi) 


1 


6 


16 


26 


31 


32 






dim^((^+) J -) 


1 


6 


16 


26 


31 


32 


32 




dini7r s (^) 


1 


6 


12 


12 


12 


12 


12 




dim^ s ((£+) x ) 


1 


6 


12 


12 


12 


12 


12 


12 



Table 9. Dimension table for tt s (Q^) and 7r s ((^) _L ) in Example 21. 

The extracted solution points 

ui = (l,8.73e-7,2.14e-<J,2.48e-7,2.23e-£,-1.29e-<5) , 

v 2 = (0.277, 0.226, 0.162, 0.0858, 0.0115, -0.124) , 

v 3 = (0.136, 0.0428, 0.0417, 0.0404, 0.0964, 0.211) , 

V4 = (0.462, 0.309, 0.0553, -0.102, -0.0844, 0.0917) , 

v 5 = (0.441, 0.151, 0.0225, 0.219, 0.0935, -0.207) , 

v 6 = (0.239, 0.0608, -0.0622, -0.0233, 0.186, 0.219) , 

v 7 = (0.753, 0.0532, 0.191, -0.114, -0.146, 0.139) , 

v s = (0.726, -0.0503, 0.122, 0.164, 0.109, -0.208) , 

v Q = (0.409, -0.0732, 0.0657, -0.127, 0.252, 0.178) , 
uio = (0.292, -0.101, 0.181, -0.0591, 0.193, 0.141) , 
vn = (0.590, 0.0422, 0.327, -0.0642, -0.0874, -0.0132) , 
v 12 = (0.68, 0.266, -0.154, 0.0323, 0.0897, -0.0735) , 

were extracted at order (t, s) — (6, 3), when conditions (3a) -(36) were first satisfied. The maximum 
evaluation error was found to be max^ e(vi) < 2.4e-^ and the commutativity error c{X) < 6.2e-6. 
Again the algorithm returns the ideal J = \fl as dim 773 (Qq) = \Vr(I)\ = 12. In this example the 
moment-matrix method [12] also extracts the 12 real solutions at order (t,s) = (6,3). 
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7. Conclusion 

This work was motivated by the great success of numerical-algebraic methods in recent years. 
Incorporating features specific to real root finding into efficient symbolic-numeric methods may 
lead to more efficient algorithms for numerically computing all real roots of a given system of 
polynomials. The contribution of this paper is a first attempt in this direction as it implements 
real-algebraic features into the existing symbolic-numeric algorithm described in [32]. Concretely, 
the resulting algorithm uses scmidcfmite programming techniques in addition to standard numerical 
linear algebra techniques. It is not only applicable to zero-dimensional ideals, but to all problems 
for which the real variety is finite. An extension to zero-dimensional basic semi-algebraic subsets 
is also possible, along the same lines as in [12]. 

The new approach relies on a dual space characterization of (an approximation of) the real 
radical ideal, obtained by combining ideas of [12] and [32], but the new prolongation-projection 
algorithm may terminate earlier than the moment-matrix method of [12]. Although preliminary 
computational results are encouraging, whether the characterization at hand can lead to a new 
treatment of real-algebraic problems is still to be demonstrated on a larger sample of problems. 
An important computational issue is how to efficiently solve the underlying semidcfinitc program 
for large problems involving high degree polynomials with many variables. Exploiting sparsity in 
order to decrease the size of the semidcfinitc program is a promising direction and the work of 
Kojima et al. [11] and Lasserre [14] is a first important step in this direction. Strategics similar to 
those used in Grobner/border basis computations can be employed to further increase efficiency of 
the proposed method, particularly in view of the linear algebra steps involved, e.g. the dimension 
tests. 
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